Optimum multiple target detection and resolution

ABSTRACT

A method or technique of resolving multiple energy sources from signals  oined by a spatial array of sensors in which portions of two prior art techniques are combined to yield an improved hybrid technique. The method involves sampling and digitizing the output of a sensor array, forming and storing a co-variance matrix from such digitized samples, determining the eigenvectors and corresponding eigenvalues relating to said matrix, performing a directional eigenanalysis on the data of the matrix and determining the minima of the resulting curves, then using a minimum variance method to calculate from said matrix data the total power at each steering direction corresponding to each of said minima, then thresholding the calculated total powers to eliminate false alarms and yield the direction of each power source.

The invention described herein may be manufactured, used and licensed by or for the Government for governmental purposes without the payment to us of any royalties thereon.

BACKGROUND OF THE INVENTION

The invention relates to the detection and resolution of multiple targets (or energy sources) from signals obtained by a spatial array of sensors. The invention has wide application in such diverse fields as passive sonar, in which case the sensors would be acoustic and the targets might be hostile submarines, earthquake and nuclear weapon detonation detection systems, in which case the sensors would be seismic and the targets the earthquake or explosion epicenters, astronomical interferometry wherein the sensors would be radiotelescopes and the targets may be distant galaxies or quasars, and phased-array radars in which case the sensors would be the array antennae. More particularly, the invention is intended to be used for precision target detection of enemy targets to enable a strike on said enemies to be commanded electronically by a battlefield command control center.

When a signal is known to consist of pure sinusoids in white noise, an appropriate procedure for determining the unknown frequencies and powers is the Pisarenko spectral-decomposition described in his paper entitled "The Retrieval of Harmonics from a Co-Variance Function", Geophysical Journal of the Royal Astronomical Society, Vol. 33, pp 347-366, 1973. All that is needed is a finite segment of the discrete covariance function whose length is at least one more than the number of sinusoids to be determined. The method is based on eigenanalysis of a Toeplitz matrix produced from the covariance function.

When Pisarenko extended this method to signals where the amplitudes and the frequencies have small random perturbations around central values he showed that linear approximation is justified when the number of sinusoids is one less than the number of eigenvectors. However, when a larger segment of the covariance function was used, leading to a larger number of eigenvectors, the statistical analysis was much more difficult. These extra eigenvectors also produced false sinusoids which appeared at fictitious frequencies with fictitious amplitudes. As a result he recommended disregarding the extra covariance data.

Since then the extension of the Pisarenko method to the problem of extracting multiple-target information from array data has led to many different techniques, e.g., the Eigen-Vector Method EVM, of Johnson and Degraff, see IEEE trans. on ASSP, Vol. 30, pp 638-647, August 1982 and the Multiple-Signal-Classification scheme (MUSIC) of Schmidt, see Proceedings RADC Sprectral Estimation Workshop, Rome, N.Y., pp. 243-258, October 1979. Between the desire to use all the available information and the problem of coping with false alarms, investigators adopted the concept of artificially dividing the eigenspace of the signal covariance matrix into two subspaces, (1) the source space, consisting of eigenvectors with large eigenvalues, and (2) noise space, consisting of the remaining eigenvectors. In these approaches the distinction between the two spaces requires a subjective judgment (or guess) on the dividing line between the sets of eigenvalues. In any case, it must be recognized that any such distinction is artificial since the noise occupies all of the eigenspace, and each eigenvector may have contributions from any source.

In these prior art procedures a metric is then formed from the so-called noise space to measure its deviation from orthogonality to the true source space. The inverse of this metric is then used to indicate the directions of arrival and the powers of the various sources. The effect of using several eigenvectors in the metric tends to decrease the number of false alarms and to produce an average location for each source. Each investigator pursuing this approach used some kind of metric, but none of them proved that the resulting procedures led to "optimum" target resolution in any sense. The present invention comprises a different approach which while based on the Pisarenko method, differs from the earlier ad hoc approaches and yields demonstrably better multi-target resolution.

SUMMARY OF THE INVENTION

The method of the present invention is a hybrid one comprising portions of two prior art techniques which when combined in the manner taught herein produces results which are superior to any prior art technique.

The concept of the present novel method of resolving multiple energy sources from received sensor array data comprises the full exploitation of the orthogonality properties of the eigenvector system, which is obtained by finding the lowest or minima of the projections of the energy sources plus noise on the eigensystem as a steering vector is scanned over the field of view. This portion of the novel method indicates the location of both true energy sources and false alarms. In accordance with the invention, the false alarms can be easily eliminated by thresholding the power at each of the candidate source locations using a known minimum variance method, such as the maximum likelihood method (MLM) of Capon, cited below, herein. The present method provides optimum source power as well as optimum source location(or resolution).

The method of the present invention comprises the steps of: forming an estimate of the covariance matrix from digitized samples (or snapshots) of the outputs of an array of sensors, performing an epgenanalysis of this matrix, determining the amount of projected energy on each eigenvector through a steering vector as that vector is steered in different directions from which target signals may come, determining the minima of said projected energy to form a set of candidate target directions which comprise true target sources as well as false alarms, then determining the total powers received from each of said candidate directions using a minimum variance method, and finally then thresholding the said power levels to eliminate the false alarms.

Another object of the invention is to provide increased resolution in the detection of targets received by a spaced array of sensors by utilizing steps from two prior art methods or techniques to form a new and improved method or technique, and wherein the novel technique of the invention does not involve dividing the signal space into noise space and source space, but involves performing an eigenanalysis on the energy contributions from all sources as well as noise associated with its eigenvector.

These and other advantages and objects of the invention will become apparent from the following detailed description and the drawings.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 is a block diagram of circuitry which may be used to carry out the method of the invention.

DETAILED DESCRIPTION OF A PREFERRED EMBODIMENT

Prior to 1969 it was commonly thought that two targets could not be resolved if the angular separation thereof were less than the beamwidth of the array of sensors. The technique used was the classical Beamformer method.

Conventional linear processing, or Beamforming, uses the response function; ##EQU1## wherein W is the steering vector corresponding to the direction θ, R is the covariance matrix of the M-element sensor array. The right hand expression represents R written in its eigenvalue expansion wherein the λ's are the eigenvalues and the V's are the eigenvectors. The superscript t indicates a conjugate transpose. Although this conventional estimator is asymtotically efficient for a single, plane wave target in white noise, its resolving power for multiple targets is limited to the received signal wavelength divided by the array maximum spatial dimension.

In 1969 Capon published his minimum variance or maximum liklihood method (MLM) in the Proceedings of the IEEE, Vol. 57, pp 1408-1418, August '69. This method uses the response; ##EQU2## This is a non-linear estimator which permits better resolution of multiple targets than that permitted by the linear response function of Equation (1).

The Johnson and De Graff (EVM) technique uses the following response; ##EQU3## Wherein k=M minus the estimated number of targets.

The MUSIC technique of Schmidt uses; ##EQU4## Wherein k has the same meaning as in Equation (3).

The aforementioned techniques, EVM and MUSIC, both are based on artificially dividing the signal space into source space and noise space wherein the noise space yields an average location of the sources if the signal to noise ratio is high enough and/or the source separation is wide enough. As a result these techniques are inadequate for many practical cases of interest.

For any sensor array the measured data can be represented as the column vector;

    X(t)=[X.sub.i (t), X.sub.2 (t), . . . X.sub.k (t)].sup.t   (5)

wherein t is time and t denotes the conjugate transpose. The cross correlation matrix of these data is;

    R=E[X·X.sup.t ]                                   (6)

wheren E here is the expectation. The eigenvectors, V, of R and their corresponding eigenvalues, λ, are defined by;

    RV.sub.k =λ.sub.k V.sub.k                           (7)

Eigenanalysis is a tool for giving us an orthogonal set of coordinates wherein the variance of the data is lowest on one axis and highest on another. These axes are the principal components.

If we now pre-multiply by V^(t) and take V^(t) V=1 we find;

    V.sub.k.sup.t RV.sub.k =E[V.sub.k.sup.t X·X.sup.t V.sub.k ]=λ.sub.k                                          (8)

which means that the variance of the projected data on a given eigenvector equals its eigenvalue. These eigenvalues can represent energy or power depending on the inuts of the data. The total energy in the system is; ##EQU5## That component of the total energy which can be projected onto any given eigenvector is its particular eigenvalue. On the other hand, the energy projected on an eigenvector V through the steering vector W is;

    φ.sub.k (θ)=λ.sub.k (W.sup.t V.sub.k V.sub.k.sup.t W) (10)

As W is steered in different directions this projection will pass through various peaks and valleys, or maxima and minima. Peaks indicate that V is near that steering direction containing a high level of received energy from a unique (single) source or possibly from an unresolved set of sources. Valleys or minima on the other hand indicate that V is nearly orthogonal to that direction of a source represented by the steering vector, or to some interference between sources.

How does one use these properties to estimate the target directions. The first guess might be to use the locations of the peaks or maxima of Equation (10), however that approach gives only slight improvement over the classical Beamformer, which looks for the maxima of the quantity; ##EQU6##

An alternate approach, which is the subject of this invention, involves performing an eigenanalysis on a covariance matrix derived from sampled and digitized snapshots of the outputs of the sensors of the array to obtain the function of Equation (10) which represents energy projected on an eigenvector through the steering vector, and hence will yield a curve of energy (or power) vs. steering angle or azimuth, the angle being related to the array geometry. The minima of this curve are selected and, as stated above, these minima represent candidate target directions which may include unresolved multiple targets and/or false alarms. Thresholding of the energy (or power) of these candidate target directions is then performed by means of a known minimum variance method, such as Capon's maximum likelihood method. This results in elimination of the false alarms.

The block diagram of FIG. 1 is acelled by load changes in another part of the system which are nearly identical in magnitude but opposite in direction. For example, the load resulting from the movement of direct sunlight around a building in the course of a day can result in changes in the load applied to various parts of a building while the total load may remain substantially the same over the full day.

Conventionally, control system response times are measured in very small fractions of a second, frequently in milli or micro seconds. However, air-conditioning load changes take place in much longer time periods as for example, periods as long as several minutes or more. When these changes take place, the control system must be sensitive enough to sense them and to react accordingly. I have recognized these characteristics and for this reason I provide a control system in which the control signal is over-dampened, thereby to increase the controller sensitivity. Over-damping is achieved by the filter 44 previously described so as to generate a signal to the controller 46 which follows a curve like curve 76 illustrated in FIG. 5 of the drawings.

The slow system response which I am able to generate is important for other reasons. Terminal boxes connected to the duct system, are usually under the control of a room thermostat. Such a thermostat will allow a typical box to pass just the correct amount of air to the space to satisfy the load. The correct amount of air depends on the pressure upstream of the box remaining constant. If the fan capacity control system which tries to keep this pressure constant is "live" and allowed to cycle, causing the duct pressure to vary, the box will allow more or less air than desirable to pass to the space causing temperature fluctuations. These fluctuations would be sensed by the thermostat which would try to correct them. The net result of this scenario is that the two control systems fight each other and the whole air-conditioning system would become unstable. The over-damping of the fan capacity controls helps to prevent this from occurring and greatly assists in providing stable controls.

As previously indicated in FIG. 2 of the drawings, the sensed variable which is detected by the detector 26 is the fan pressure. Generally, this pressure is not steady enough for control purposes and the filter 44 which has a long time constant serves to filter out transients and to allow only actual load changes heavily filtered to pass to the controller 46. The time constant is adjustable with the adjustment providing a range of 0 to 200 seconds. The time constant may therefore be adjusted on-site to be a good compromise between the desired system response time and control system stability.

As previously indicated, integrators 62 or 64 may be substituted for or used in conjunction with the filter 44. In control systems requiring fast response times, cycling occurs at a specific frequency. If the control system is set for the desired sensitivity and impermissable and excessive cycling occurs, substituting an integrator for a filter may help to overcome this difficulty. The integrator must be inserted into the circuit as illustrated in FIG. 3 of the drawings between the signal conditioner 42 and the controller 46. As the transmitter senses the system cycling, the integrator will act on the output. As long as the cycling is of a set frequency, the integrator will integrate out the fluctuations. The integral of a symetrical periodic function is zero. Therefore, the integrator output will consist of the steady state value only as sensed by the transmitter. It is important that the frequency of the cycling be measured on the job site and adjusted into the integrator. An alternate position of the integrator is in the feedback loop. This device may be used with or without filter 44 and integrators 62 and 64.

Controller 30 may also be constructed with a micro-processor replacing some or all signal processing, conditioning, filtering, integrating, controlling, feedback, set point and detecting functions by means of software algorithms contained in the micro-processor's memory.

In air delivery systems having well designed air entry and exit configurations for the fan, the accuracy of my method and control system is derived from the accuracy of the fan curves which are provided by the manufacturers of the various fans. The manufacturers take great pains to assure the accuracy of their fan curves. It will be apparent that the manufacturer's fan curves may have to be modified for a particular installation if the inlet and outlet configuration of the fan are significantly different from those utilized by the fan manufacturer in preparing the ideal fan curves.

While some inaccuracy will inevitably be present in my system this will primarily result from loading of the air or fluid filters which are employed in my system. However, if the pressure drop across the filters due to clogging increases by 1 inch W.G. in a fan system selected for 4 inches W.G., the error only increase by 3%. Consequently the accuracy of my system is improved by a significant factor. Furthermore, because all signalling processes are done electronically, drifting is essentially eliminated and recalibration is essentially no longer required.

Various modifications of the present invention will be apparent to those skilled in the art. It will be apparent that the air conditioning control system of the present invention may be used to advantage in a single storey building in which case it serves to overcome the difficulties experienced with lateral rather than vertical air disturbances. 

I claim:
 1. In a fluid flow system which includes a conduit for conveying the fluid from a fluid delivery device which provides a variable volume fluid discharge by varying the speed of operation of the fluid delivery device, and wherein it is necessary to generate a rate of flow signal which is proportional to the volume of fluid delivered by the fluid delivery device in use, the improved method of generating said rate of flow signal comprising the steps of;(i) generating a first input signal which is proportional to the speed of the fluid delivery device, (ii) inputting the first input signal to a curve fitter which programmed with the combined fluid delivery device/conduit system characteristic curve according to the formula;

    C.sub.1 q.sup.2 -KS.sup.2 =0

where;C₁ is a constant (turbulent flow coefficient) q is the rate of flow in the system, K is constant =P₂ /S₂ ² where;P₂ is the maximum pressure generated by the fluid delivery device when it operates at design conditions, S is the speed of the fluid delivery device, S₂ is the speed of the fluid delivery device when operating at maximum design conditions and to thereby generate said rate of flow signal, as an output signal from said curve fitter, which is proportional to fluid flow delivered by the fluid delivery device.
 2. The method of claim 1 wherein said fluid delivery device is a variable speed fan.
 3. In a fluid flow system which includes a conduit for conveying the fluid from a fluid delivery device which provides a variable volume fluid discharge by varying the speed of operation of the fluid delivery device, and wherein it is necessary to generate a rate of flow signal which is proportional to the volume of fluid delivered by the fluid delivery device in use, the improved method of generating said rate of flow signal comprising the steps of;(i) generating a first input signal which is proportional to the speed of the fluid delivery device, (ii) generating a second signal which is proportional to the pressure in the conduit at zero flow, and (iii) inputting said first and second input signals to the first and second inputs respectively of a dual input curve fitter which is programmed with the combined fluid delivery device/conduit system characteristic curve according to the formula;

    C.sub.1 q.sup.2 +C.sub.3 q.sup.0 -KS.sup.2 =0

where; C₁ is a constant (turbulent flow coefficient) C₃ is a constant which is independent of flow and is equal to the pressure generated by the fluid device at zero flow.q is the rate of flow in the system, K is a constant=P₂ /S₂ ² where;P₂ is the maximum pressure generated by the fluid delivery device when it operates at design conditions, S is the speed of the fluid delivery device, S₂ is the speed of the fluid delivery device when operating at maximum design design conditions and to thereby generate said rate of flow signal, as an output signal from said curve fitter, which is proportional to fluid flow delivered by the fluid delivery device.
 4. The method of claim 3 wherein said fluid delivery device is a variable speed fan.
 5. In a fluid flow system which includes a conduit for conveying the fluid from a fluid delivery device which provides a variable fluid discharge by varying the speed of operation of the fluid delivery device, and wherein it is necessary to generate a rate of flow signal which is proportional to the volume of fluid delivered by the fluid delivery device in use, the improved method of generating said rate of flow signal comprising the steps of;(i) generating a first input signal which is proportional to the speed of the fluid delivery device, (ii) generating a second signal which is proportional to the pressure in the conduit at zero flow, and (iii) inputting said first and second input signals to the first and second inputs respectively of a dual input curve fitter which is programmed with the combined fluid delivery device/conduit system characteristic curve according to the formula;

    C.sub.1 q.sup.2 +C.sub.2 q.sup.1 +C.sub.3 q.sup.0 -KS.sup.2 =0

where;C₁ is a constant (turbulent flow coefficient) C₂ is a constant (laminar flow coefficient) C₃ is a constant which is independent of flow and is equal to the pressure generated by the fluid device at zero flow. q is the rate in the system, K is a constant=P₂ /S₂ ² where;P₂ is the maximum pressure generated by the fluid delivery device when it operates at design conditions, S is speed of the fluid delivery device, S₂ is the speed of the fluid delivery device when operating at maximum design conditions and to thereby generate said rate of flow signal, as an output signal from said curve fitter, which is proportional to fluid flow delivered by the fluid delivery device.
 6. The method of claim 5 wherein said fluid delivery device is a variable speed fan. 